I wanted to look at absorption features from mineral spectral signatures to determine type of hydrothermal alteration and what minerals you can expect to find just through remote data.
Reading in packages and all of the absorption feature plots I saved
library(terra)
library(raster)
library(future)
library(future.apply)
library(doParallel)
library(zoom)
library(tidyverse)
library(janitor)
library(grid)
library(gridExtra)
library(circular)
source("../mineral_research/functions.R")
area_750_1250 <- (rast("../AVIRIS/2023/AV320231005t185153/area_750_1250.tif"))
area_1300_1600 <- (rast("../AVIRIS/2023/AV320231005t185153/area_1300_1600.tif"))
area_1600_2200 <- (rast("../AVIRIS/2023/AV320231005t185153/area_1600_2200.tif"))
area_1750_2100 <- (rast("../AVIRIS/2023/AV320231005t185153/area_1750_2100.tif"))
area_2050_2250 <- (rast("../AVIRIS/2023/AV320231005t185153/area_2050_2250.tif"))
area_2150_2250 <- (rast("../AVIRIS/2023/AV320231005t185153/area_2150_2250.tif"))
area_2275_2375 <- (rast("../AVIRIS/2023/AV320231005t185153/area_2275_2375.tif"))
area_2250_2450 <- (rast("../AVIRIS/2023/AV320231005t185153/area_2250_2450.tif"))
area_2200_2400 <- (rast("../AVIRIS/2023/AV320231005t185153/area_2200_2400.tif"))# Define your desired extent
desired_extent <- ext(402000, 404500, 4246000, 4248000)
#reading in 2023 AVIRIS
cropped_153 <- (rast("../AVIRIS/2023/AV320231005t185153/L2A_OE/cropped_1.tif"))
#cropping
cropped_153 <- crop(cropped_153,desired_extent)
#reading in saved scores from PCA
scores <- rast("../AVIRIS/2023/AV320231005t185153/scores.tif")
#reading in 2019 AVIRIS
aviris_2019 <- rast("../AVIRIS/2019/rectified.tif")
#cropping 2019 AVIRIS
aviris_2019 <- crop(aviris_2019,desired_extent)
#performing PCA
PCA2 <- prcomp(aviris_2019)
#turning into rasater of PCA
scores2 <- predict(aviris_2019,PCA2,index=1:5)
# Get the CRS of scores2
crs_scores2 <- crs(scores2)
#reading in geologic map
malmsten_cropped <- rast("../Quads/Malmsten_peak_quad/cropped_malmsten.tif")
#reading in shape file of geologic units
shape <- vect("../Quads/Malmsten_peak_quad/units_malmsten.shp")
#projecting shape units to the same reference system
shape <- project(shape,crs_scores2)
#cropping shapefile
shape_cropped <- crop(shape,desired_extent)
# Specify the directory containing the .tif files
directory <- "../landsat/LC09_L2SP_038033_20231019_20231020_02_T1/LC09_L2SP_038033_20231019_20231020_02_T1/"
# Get a list of all .tif files in the directory
tif_files <- list.files(directory, pattern = "\\.TIF$", full.names = TRUE)
tif_files <- rast(tif_files)
tif_files <- crop(tif_files,desired_extent)
PCA3 <- prcomp(tif_files)
scores3 <- predict(tif_files,PCA3)Plotting the PCA for the AVIRIS-3 image
# Define custom colors for each unit
unit_colors <- c(Qa = "cornsilk1", Ql = "burlywood1", Qp = "lightgreen", Tdk = "hotpink", Tlv = "red4", Tn = "red", Tda = "purple",Tla = "darkcyan")
# Plot the RGB image
plotRGB(scores,1,2,3, axes=TRUE, stretch="lin", main="2023 AVIRIS (3m resolution) PCA RGB")Plotting the PCA for the AVIRIS-NG image
Plotting the digitized geologic map
plot(malmsten_cropped, main = "Geologic map")
plot(shape_cropped, col = unit_colors[shape_cropped$unit],add=T)
legend("topleft", legend = names(unit_colors), fill = unit_colors, cex = 0.7,inset = c(0, 0.3))Reading in spectral data from the USGS and cleaning it
#reading in spectral mineral signatures
list <-
list.files(path ="../usgs_spectral_library/usgs_splib07 (1)/ASCIIdata/ASCIIdata_splib07b_cvAVIRISc2014/ChapterM_Minerals/",
full.names = TRUE)
#list applying as csv and into a data frame
minerals <- list %>%
lapply(read.csv) %>%
as.data.frame() %>%
clean_names() %>%
mutate(Bands = row_number())
#removing everything before first underscore
colnames(minerals) = sub(".*?_", "", colnames(minerals))
#removing everything before first underscore
colnames(minerals) = sub(".*?_", "", colnames(minerals))
#removing everything before first underscore
colnames(minerals) = sub(".*?_", "", colnames(minerals))
#keeping only after first two underscores
colnames(minerals) <- str_extract(colnames(minerals), "[^_]*_[^_]*[^_]*")
#renaming blank column to band
colnames(minerals)[1277] = "Band"
# Generate the vector
vector <- seq(from = 400, by = 9.375, length.out = 224)
#adding wavelength as a column
minerals <- minerals %>%
mutate(wavelength = vector)Plotting a bunch of mineral spectral signatures from different hydorthermal alteration types
# Identify columns with negative values
neg_cols <- sapply(minerals, function(x) any(x < 0))
# Subset the data, excluding columns with negative values
data_subset <- minerals[, !neg_cols]
# Define custom breaks for the x-axis
custom_breaks <- seq(400, 2500, by = 100)
#chlorite spectral signature
data_subset %>%
ggplot(aes(x=wavelength)) +
geom_line(color="blue",aes(y=data_subset$'805_chlorite'),size=1) +
geom_line(color="blue",aes(y=data_subset$'1140_epidote'),size=1)+
geom_line(color="green",aes(y=data_subset$'2206_muscovite'),size=1)+
geom_line(color="green",aes(y=data_subset$'1768_kaolinite'),size=1)+
geom_line(color="red",aes(y=data_subset$'1412_halloysite'),size=1)+
geom_line(color="red",aes(y=data_subset$'1031_dickite'),size=1)+
geom_line(color="orange",aes(y=data_subset$'2788_quartz'),size=1)+
geom_line(color="orange",aes(y=data_subset$'742_chalcedony'),size=1)+
geom_line(color="purple",aes(y=data_subset$'2736_pyrophyllite'),size=1)+
geom_line(color="purple",aes(y=data_subset$'1020_diaspore'),size=1)+
geom_line(color="yellow",aes(y=data_subset$'581_biotite'),size=1)+
geom_line(color="black",aes(y=data_subset$'90_albite'),size=1)+
geom_line(color="black",aes(y=data_subset$'2341_nepheline'),size=1)+
ylim(0,1) +
xlim(420,2500) +
labs(y="Reflectance",
x="Wavelength (nm)") +
scale_x_continuous(breaks = custom_breaks) +
theme_bw()#Albite spectral signature
data_subset %>%
select(matches("albite|wavelength")) %>%
pivot_longer(cols = ends_with("albite")) %>%
filter(name %in% unique(head(name, 10))) %>% # Keep data for the first 8 unique chlorites
ggplot(aes(x = wavelength, y = value, col = name)) +
geom_line() +
labs(x = "Wavelength (nm)",
y = "Reflectance",
col = NULL) +
theme_bw()This is how I determined the mean area and standard deviation for all the minerals in the wavelength ranges.
#Pivoting longer
data_subset_long <- pivot_longer(data_subset,
cols = -c(Band, wavelength),
names_to = "Sample_Mineral",
values_to = "Reflectance")
# Extract sample and mineral names
data_subset_long <- separate(data_subset_long, Sample_Mineral, into = c("Sample", "Mineral"), sep = "_")
# Convert wavelength column to numeric
data_subset_long$wavelength <- as.numeric(data_subset_long$wavelength)
# Renaming the 'wavelength' column to 'Wavelength'
data_subset_long <- rename(data_subset_long, Wavelength = wavelength)
visualize_area_df_mineral(data_subset_long, mineral_name = "halloysite",wavelength1 = 1300,wavelength2 = 1600)## $plots
## $plots$`1412`
##
## $plots$`1418`
##
## $plots$`1423`
##
## $plots$`1430`
##
## $plots$`1428`
##
##
## $areas_df
## Sample Area
## 1 1412 13.45025
## 2 1418 17.24728
## 3 1423 19.85488
## 4 1430 20.43367
## 5 1428 13.86969
# Call the function with your data
results <- visualize_area_all_minerals(data_subset_long, wavelength1 = 750, wavelength2 = 1250)
#arranging in decending order
results <- arrange(results,desc(Area))
# Assuming your data frame is named 'df'
mean_area_by_mineral <- results %>%
group_by(Mineral) %>%
summarize(mean_area = mean(Area),
sd_area = sd(Area))
# Print the result
head(arrange(mean_area_by_mineral,desc(mean_area)),n =5)## # A tibble: 5 × 3
## Mineral mean_area sd_area
## <chr> <dbl> <dbl>
## 1 bronzite 36.8 8.34
## 2 jarosite 31.4 7.86
## 3 hypersthene 30.7 11.4
## 4 rhodonite 27.3 NA
## 5 lepidocrocite 26.8 NA
Reading in dem and creating slope and aspect maps
#reading in dem over the area
dem <- rast("../DEMs/dry_creek/output_USGS1m.tif")
#cropping by extent
dem <- crop(dem,desired_extent)
#creating slope and aspect map
slope <- terrain(dem, v = "slope", unit = "degrees")
aspect <- terrain(dem, v = "aspect", unit = "degrees")Sobel Filtering
# Calculate gradients (Sobel filters)
gx <- focal(dem, matrix(c(-1, 0, 1, -2, 0, 2, -1, 0, 1), nrow = 3), fun = sum)
gy <- focal(dem, matrix(c(1, 2, 1, 0, 0, 0, -1, -2, -1), nrow = 3), fun = sum)
edges = sqrt(gx^2 + gy^2)
# Apply threshold to gradients to identify significant features
threshold <- 20 # Adjust threshold as needed
gx[edges <= threshold] <- NA
gy[edges <= threshold] <- NA
# Convert edge data to slope angles (assuming they represent the orientations)
lineament_angles <- atan2(na.omit(values(gy)), na.omit(values(gx))) * 180 / pi
# Adjust angles to be in the range of 0 to 360 degrees
lineament_angles <- (lineament_angles + 360) %% 360
# Convert angles to radians for circular statistics
lineament_angles_rad <- lineament_angles * pi / 180
# Restrict angles to the range 0-180 and 180-360
lineament_angles[lineament_angles > 180] <- lineament_angles[lineament_angles > 180] - 180
# Create circular data object
angles_circ <- circular(lineament_angles, units = "degrees")
# Define the breaks for grouping the angles
breaks <- seq(0, 360, by = 360 / 24)
# Use the cut function to group the angles into bins
angle_bins <- cut(lineament_angles, breaks = breaks, include.lowest = TRUE)
# Calculate the frequency of each angle bin
angle_counts <- table(angle_bins)
# Calculate the maximum count to use as the reference for scaling
max_count <- max(angle_counts)
# Scale the lengths of the bars based on the frequency of each angle bin
scaled_lengths <- angle_counts / max_count
plot(shade(slope*pi/180, aspect*pi/180, angle = 10, direction = 270), col = grey(0:100/100), legend = F)
plot(edges > threshold, col = c("transparent","red"), add = T)# Plot the rose diagram with adjusted lengths
rose.diag(angles_circ, bins = 24, prop = 2, col = "red")
rose.diag(angles_circ, bins = 24, col = "red",prop = 2, zero =pi, add = T, axes = F)Making a map from all the minerals that showed up as having the most absorption area for the wavelength ranges with one standard deviation.
# Define colors and associated minerals
colors <- c("red", "orange", "green", "blue", "purple")
minerals <- list(
c("Analcime", "Halloysite", "Saponite", "Heulandite"),
c("Dickite", "Nacrite", "Paragonite", "Kaolinite"),
c("Bronzite", "Jarosite", "Hypersthene", "Rhodonite"),
c("Analcime", "Saponite", "Laumontite", "Clinoptilolite"),
c("Zunyite", "Dickite", "Alunite", "Pyrophyllite")
)
plot(shade(slope*pi/180, aspect*pi/180, angle = 10, direction = 270), col = grey(0:100/100), legend = F)
# Add other plots with 'add = TRUE' to overlay them
plot(area_1300_1600, breaks = c(17-3.26,17,17+3.26), col = c("red","red"), legend = FALSE, add = T)
plot(area_2150_2250, breaks = c(5.47-1.06,5.47,5.47+1.06), col = c("orange","orange"), add = TRUE, legend = FALSE)
plot(area_750_1250, breaks = c(37-8.34,37,37+8.34), col = c("green","green"), add = TRUE, legend = FALSE)
plot(area_1750_2100, breaks = c(22.3-1.45,22.3,22.3+1.45), col = c("blue","blue"), add = TRUE, legend = FALSE)
plot(area_2050_2250, breaks = c(8.64-2.68,8.64,8.64+2.68), col = c("purple","purple"), add = TRUE, legend = FALSE)
# Plot colored rectangles with mineral labels
legend_x <- 402000
legend_y <- 4247800
legend_rect_size <- 80
legend_text_x <- legend_x + legend_rect_size + 50
for (i in seq_along(minerals)) {
legend_labels <- minerals[[i]]
legend_color <- colors[i]
legend_rect_y <- legend_y - i * legend_rect_size
rect(legend_x, legend_rect_y, legend_x + legend_rect_size, legend_rect_y - legend_rect_size, col = legend_color)
text(legend_text_x, legend_rect_y - legend_rect_size / 2, paste(legend_labels, collapse = ", "), adj = 0)
}Reading in site collection location data and the spectrometer images.
#reading in sample location sites
locat <- read.csv("./Locations_minerals.csv")
locat_v <- vect(locat, geom = c("Longitude", "Latitude"), crs = "EPSG:4326")
plot(shade(slope*pi/180, aspect*pi/180, angle = 10, direction = 270), col = grey(0:100/100), legend = F, mar = 3)
locat_v_pro <- project(locat_v,slope)
plot(locat_v_pro, "Title", add = T, col = c("red","orange","green","blue","purple","hotpink"))Sample 1
Sample 2
Sample
3
Sample 4
Sample 5
Sample 6
Here is the scolecite mineral signature from the USGS, it looks very similar to the sample 5!
#Scolecite spectral signature
data_subset %>%
select(matches("scolecite|wavelength")) %>%
pivot_longer(cols = ends_with("scolecite")) %>%
ggplot(aes(x=wavelength,y=value,col = name)) +
geom_line() +
labs(x="Wavelength (nm)",
y="Reflectance") +
theme_bw()This project concluded that there is zeolitic alteration (a low grade type of hydrothermal alteration) in this area, due to the moderate welding of the andesite and the presence of Scoelecite.